# -*- coding: utf-8 -*-

"""
    ** NOTE: SOLUTION IS NOT CORRECT **

    http://projecteuler.net/problem=394
    
    PROBLEM
    
    Jeff eats a pie in an unusual way.
    
    The pie is circular. He starts with slicing an initial cut in the pie
    along a radius.
    
    While there is at least a given fraction F of pie left, he performs the
    following procedure:

    -   He makes two slices from the pie centre to any point of what is
        remaining of the pie border, any point on the remaining pie border
        equally likely. This will divide the remaining pie into three pieces.

    -   Going counterclockwise from the initial cut, he takes the first two pie
        pieces and eats them.
        
    When less than a fraction F of pie remains, he does not repeat this
    procedure. Instead, he eats all of the remaining pie.

    For x ≥ 1, let E(x) be the expected number of times Jeff repeats the
    procedure above with F = 1/x.

    It can be verified that E(1) = 1, E(2) ≈ 1.2676536759, and
    E(7.5) ≈ 2.1215732071.

    Find E(40) rounded to 10 decimal places behind the decimal point. 

    
    NOTES
    Wasn't able to crack this. Current solution:
    
    x   E(x)            Expected
    2   1.2656288147    1.2676536759
    7.5 2.25139622703   2.1215732071
    
    
    REFERENCES
    http://math.stackexchange.com/q/72589
    http://en.wikipedia.org/wiki/Expected_value
    
    
    PERFORMANCE

    
"""
#
# Import
#
import time
import math, random
from numpy import *
from decimal import Decimal


#
# Globals / Constants
#
F1 = 1.0
D1 = Decimal(1.0)
MIN_SLICE = 10**-10

def timeit(f):
    def timer():
        t0 = time.time()
        returned = f()
        print "time %s: %.6f s" % (f, time.time() - t0)
        return returned
    return timer

def assert_match(value, expected, message=""):
    if not message:
        message = "error: %.9f %%" % ( error_margin(value, expected) * 100 )
    assert value == expected, "value %s != expected %s (%s)" % (
        value, expected, message)
    
def error_margin(value, expected):
    return (float(value) - expected) / expected


#
# Test Case / Solution
#
@timeit    
def test_case():
    cases = [
        ([.24, .43, .68, .89], 0.98475136),
        ([.75, .96], 0.99),
        ([.75], 0.75),
    ]
    #for l, e in cases:
    #    print pN(l), e
    
    cases = [
        (2.0, 1.2676536759),
        (7.5, 2.1215732071),
    ]
    
    for x, expected_slices in cases:
        slices = E(x)
        print x, slices, expected_slices, "\n\n"

    assert_match(slices, expected_slices)
    
    print "test case passed!"

@timeit
def solution():
    g = euler_accelerator(G(40))
    [g.next() for _ in xrange(1000000)]
    return g.next()
        


#
# Solution Code
#
def E(x, ps=None):
    if not ps:
        ps = []
        
    p = p2(1.0/x, 1.0/x) * (1.0 - sum(ps))
    ps += [p]
    x = 1.0 - sum(ps)
    print len(ps), p, x, sum(ps), ps
    
    if p < MIN_SLICE:
        return sum((n*p) for n,p in enumerate(ps,1)) / sum(ps)
    else:
        return E(x, ps)
        
def slice_it(x, p):
    edible = 1.0 / x
    poison = 1.0 - edible
    
    sliced = poison * p
    poison = poison - sliced
    
    inv_x = edible / (edible + poison)
    return 1.0 / inv_x
    

def Z(x, ps=None):
    if not ps:
        ps = []
    
    p = 1.0/x
    x = p2(p, p)
    xa = x * (1.0 - sum(ps))
    #xa = x - sum(ps)
    ps += [xa]
    print len(ps), x, ps
    x = 1.0 / p2(p, x)
    
    if 1.0 - sum(ps) < MIN_SLICE:
        return sum((n*p) for n,p in enumerate(ps,1))
    else:
        return Z(x, ps)
    

def AR(x, slices=0, psum=1.0):
    edible = 1.0 / x
    sliceable = 1.0 - edible
    slices += 1
    
    p = p_slice_is_last(sliceable, edible)
    p_adj = psum * p
    psum -= p_adj
    
    #slice_size = p * sliceable
    #sliceable -= slice_size
    #x = (edible + sliceable) / edible
    edible = p
    sliceable = 1.0 - p
    x = (edible + sliceable) / edible
    
    
    print slices, p, p_adj, psum, x, (edible, sliceable)
    
    if p_adj <= MIN_SLICE:
        return slices * p_adj
    else:
        return (slices * p_adj) + R(x, slices, psum)

def A(x):
    """adjust sliceable edible to probs"""
    pie = 1.0
    edible = pie / x
    sliceable = pie - edible
    slices_adj = 0

    p_adj = MIN_SLICE
    p_run = 0
    n = 0
    ps = []
    pas = [0]
    
    while p_adj >= MIN_SLICE:
        n += 1
        p = p_slice_is_last(sliceable, edible)
        p_adj = pie * p
        pie -= p_adj
        slice_size = p * sliceable
        slices_adj += n * p_adj
        ps += [p]
        pas += [p_adj]
        
        # p_run
        #p_run = pN(ps)
        
        edible = p
        sliceable = 1.0 - p
        
        print n, p_adj, p, slices_adj
    
    return slices_adj

def D(x):
    """A with decimals"""
    x = Decimal(x)
    pie = D1
    edible = pie / x
    sliceable = pie - edible
    slices_adj = 0

    p_adj = Decimal(MIN_SLICE)
    p_run = 0
    n = 0
    ps = []
    pas = [0]
    
    while p_adj >= Decimal(MIN_SLICE):
        n += 1
        p = d_slice_is_last(sliceable, edible)
        p_adj = pie * p
        pie -= p_adj
        slice_size = p * sliceable
        slices_adj += n * p_adj
        ps += [p]
        pas += [p_adj]
        
        # p_run
        #p_run = pN(ps)
        
        edible = p
        sliceable = D1 - p
        
        print n, p_adj, p, slices_adj
    
    return slices_adj


def R(x, slices=0, psum=1.0):
    edible = 1.0 / x
    sliceable = 1.0 - edible
    
    slices += 1
    p = p_slice_is_last(sliceable, edible)
    p_adj = psum * p
    psum -= p_adj
    
    slice_size = p * sliceable
    sliceable -= slice_size
    x = (edible + sliceable) / edible
    
    print slices, p, p_adj, psum, x, (edible, sliceable)
    
    if p_adj <= MIN_SLICE:
        return slices * p_adj
    else:
        return (slices * p_adj) + R(x, slices, psum)
        
 
def p_slice_is_last(sliceable, edible):
    pie = sliceable + edible
    p_last = edible / pie
    return p2(p_last, p_last)
    
def d_slice_is_last(sliceable, edible):
    pie = sliceable + edible
    p_last = edible / pie
    return pd(p_last, p_last)
    
def pd(p1, p2):
    return p1 + ((D1 - p1) * p2)

def p2(p1, p2):
    """probability of p1 or p2 happening at least once"""
    return p1 + p2 - (p1 * p2)
    
def pN(ps, first_p=None):
    if len(ps) == 1 and first_p:
        return p2(first_p, ps)
    elif len(ps) == 1:
        return ps[0]
    else:
        return p2(pN(ps[1:]), ps[0])
        

def G(x):
    slices = 0
    trials = 0
    
    while True:
        slices += jeff_eats(x)
        trials += 1
        yield float(slices) / trials
        

def jeff_eats(n):
    pie = math.pi * 2
    eat_point = pie / n
    slices = 0
        
    while pie > eat_point:
        pie = jeff_slices(pie)
        slices += 1
        
    return slices
    

def jeff_slices(pie):
    return pie - (pie * max([random.random(), random.random()]))
    
    
def euler_accelerator(g):
    s0 = g.next() # Sn-1
    s1 = g.next() # Sn
    s2 = g.next() # Sn+1
    while 1:
        try:
            yield s2 - ((s2 - s1)**2)/(s0 - 2*s1 + s2)
            s0, s1, s2 = s1, s2, g.next()
        except:
            yield g.next()
            
            
def polyfitter(n, dataset, degree=3):
    import numpy
    
    xs = len(dataset) * [n]
    ys = dataset
    co = numpy.polyfit(xs, ys, degree)
    
    eq = lambda x,a: sum(c*(x**(len(a)-i-1)) for i,c in enumerate(a))
    
    return eq(n, co)

    

#
# Additional Tests
#
def test_brute_performace():
    precision = 10**9
    limit = 100000
    t0 = time.time()

    for n in xrange(limit):
        jeff_eats(7.5)
        
    secs = time.time() - t0
        
    rate = limit / secs
    est_time = precision / rate
    print rate, est_time / 60, "mins"
    
def test_slice_limit():
    n = 100000
    print sum(max([random.random(), random.random()]) for _ in xrange(n)) / n
            

#
# Main
#
if __name__ == '__main__':
    #test_brute_performace()
    #test_slice_limit()    
    test_case()
    print solution()
    